FINAL REPORT 



Nonlinear Acoustic Wave Propagation in Atmc 


Principal Investigator: Dr.S.l.Hariharan 


Period: October 1, 1984 through October 31 1985 

Granting Institution: The University of Tennessee Space Institute 

Tullahoma, TN 37388 

Grant Number: NAG - 1 - 527 


[NASA-CR- 176642} NON LINE AR ACOUSTIC WAVE N86-22309 

PROPAGATION IN ATEOSPEEfiE Final Report, 1 
Oct. 1984 - 3 1 Cct. 1 S £ 5 (Tennessee Univ. 

Space Inst., Tullahoma.) 28 p EC A03/MF A01 Onclas 

CSCL 20A G3/71 16954 



Nonlinear Acoustic Wave Propagation in Atmosphere 


S.I.Hariharan 

Department of Mathematical Sciences 
University of Akron, Akron, OH 44325 

ABSTRACT 

In this paper we consider a model problem that simulates an 
atmospheric acoustic wave propagation situation that is nonlinear. The 
model is derived from the basic Euler equations for the atmospheric 
flow and from the regular perturbations for the acoustic part. The 
nonlinear effects are studied by obtaining two successive linear 
problems in which the second one involves the solution of the first 
problem. Well-posedness of these problems is discussed and 
approximations of the radiation boundary conditions that can be used in 
numerical simulations are presented. 


* Research for the author was supported by a grant from NASA, Grant 
No. NAG 1-527 through University of Tennessee Space Institute, 
Tullahoma, TN, and by a NASA contract no. NAS 17070 while the 
author was in residence at ICASE , NASA Langley Research Center, 
Hampton, VA. 


i. Introduction 


In this paper we are interested in a two dimensional model of 
acoustic wave propagation in the atmosphere. The propagation originates 
from a point source with a high intensity of sound. It is well known that 
acoustic wave propagation in the atmosphere is rather a complex 
phenomenon. It is influenced by atmospheric conditions such as 
pressure, density, temperature, and wind variations. To analyze the 
complete problem is a difficult task. However, numerical methods have 
proven capabilities of handling such problems, but it has not yet been 
carried out for this class of problems. During the 1960’s approximate 
analytical methods have been attempted for simplified models of the 
atmosphere. Axisymmetric three dimensional models were done by Cole 
and Greifinger [1], and discussions on the physical nature were done by 
Pierce [2]. These models essentially handle only linear wave 
propagation but allow inhomogeneities in the atmospheric conditions on 
pressure and density. Our ultimate goal is to treat the full nonlinear 
model which can incorporate all variations of atmospheric conditions. 
However, at the present time we shall be concerned with the simplified 
situations of the above model, but retaining nonlinearity. 

The goal of the paper is two- fold. First one is to define the problem 
that governs the nonlinear behavior. It turns out that the problem can be 
decomposed into two linear ones. We examine the well-posedness of 
these problems, i.e., devise the mechanism that will lead to the 
existence and uniqueness of the solution of these problems. The second 
part will contain a brief discussion on radiation or absorbing conditions 
that are suitable for numerical calculations. It turns out that 
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these are derived from the appropriate dispersion relations for the 
linear problem. The question of well-posedness of the problem also 
plays a crucial role in the numerical simulation. In this particular 
study, modelling the acoustic source is a crucial part. The source 
should be of the nature that does not violate the well-posedness of the 
problem. Here we are interested in including sources that are rather 
nonsmooth. An example is a blast wave sound such as a space shuttle 
takeoff situation. The sound sources may be pulses, i.e., " delta 
functions" in both space and time. However, other standard sinusoidal 
types of sources can be included without difficulty. It turns out that the 
well posedness of the problem yields the regularity of the source and 
thus gives a guideline to approximate a source such as the delta 
function type in a proper manner. The analysis of this class of 
problems can be treated according to the theory of Kreiss[3] for 
hyperbolic systems. Unfortunately what turns out is a characteristic 
problem. This does not conform with Kreiss’s analysis entirely. We 
take a slight deviation from his approach. As an outcome of this 
analysis, one can also derive a family of boundary conditions that can be 
used for numerical computations. The numerical results will be 
reported elsewhere. 

As we mentioned earlier the governing equations are derived from 
the Euler equations for the atmosphere. As in Cole and Greifinger [1], 
we consider an isothermal atmosphere, above a ground plane with sound 
produced by instantaneous energy release at a point on the ground. We 
will also consider cases other than that of instantaneous release rate 
such as sources of smooth sinusoidal type. 
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We shall begin with the statement of the fluid flow problem which 
governs the acoustic phenomena. The model assumes strictly a two 
dimensional field of atmosphere with a source that produces a release 
of sound energy at a point Zq above the ground. The case we treat is of 
an isothermal atmosphere with the standard model of exponentially 
varying pressure and density fields. We shall not be concerned with 
wind speed so that the atmosphere is in equilibrium. Then the 
equilibrium atmosphere is characterized by exponential distributions 
for pressure and density with a scale height h, 




( 1 . 1 ) 


where P, p, and T are sea level pressure, density, and temperature 

respectively. Also, the scale height h is given by 

★ 


h = 


RT 

T 


( 1 . 2 ) 


To nondimensionalize the problem we need 

c*= V yRT* (1.3) 


which is the isentropic sound speed and 

c g - Vgh (1.4) 
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the gravity sound speed. Then the nondimensional form of the Euler 
equations ( the equations of continuity, balance of momentum, and 
energy) is 


dp 


Si +div (pg) 

= 0 

(1.5) 
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Note that in equation (1.6) the forcing term, -i_ k (k is the unit vector 

T 

in the z direction) , arises due to the forcing term per unit mass -gk in 
the original variables which is due to gravity. In equation (1.7), f(x,z,t) 
dictates the space time dependency of the source (see figure 1) and e 
measures the energy release per unit volume. For the case of an 
instantaneous energy release, e is given by 


e 


= (y-i)Q 0 



(1.8) 


where Qq is the total energy released at time t = 0. The initial 
conditions are 

p — p — e z , q = 0 at t = 0. (1.9) 

The boundary conditions at z = 0 are 
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( 1 . 10 ) 


a z = g 

which states that the vertical component of the flow is zero at z = 0. 


source ) 


I" {0> 2 q) 


////////////////////// < ground 


2. Formulation of tiie Acoustics Problem 
The acoustic expansion is based on e << 1 and represents the flow as 
small changes superimposed on the flow of the ambient state. We note 
that the ambient velocity is zero , but pressure and density have the 
form e . z Thus the expansions are 

g = eu + e^Uj + (2.1) 

p = e z { 1 + ep + e^pj + } (2.2) 

p = e z { 1 + ecr + e^cr^ + } (2.3) 

where u = (u,w) , and u^ = (UpWj). Quantities u,u^ and w,Wj are x 
and the z components of the acoustic velocities, respectively. We 
substitute expansions (2.1) - (2.3) into equations (1.2) -(1.4), initial 
conditions (1.7) and boundary conditions (1.8) to obtain the field 
equations. 

Problem that results from order e is linear and is similar to the 
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one reported by Cole and Greifinger [1]. This is as follows: 


<7, + U + W - W = 0 
t X z 

(2.5) 

u t + !p x =0 
y 

(2.6) 

w +J_p -£0; =0 

y y 

(2.7) 

P t - y<7 t + (y-i)w = f(x,z,t). 
We rewrite (2.8) using (2.5) as 

(2.8) 

p t + ^ + yw z " w = f(x,z,t) . 

(2.8) 

Initial and boundary conditions for these perturbations are 
p = cr = u = w = 0 

(2.9) 

w = 0 at z = 0 , t > 0 . 

(2.10) 

For convenience we will call (2.5) - (2.10) problem Pj. 
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Similarly the problem that results from order e has the form 

<T l,t + u l > x + w l, Z - w l = f l 

(2.11) 

u l,t + (1/y) p 1>x = f 2 

(2.12) 

w l,t + P2,z" = f 3 

(2.13) 

Pl,t'7 <J l,t +( Tf 1)w l = f 4- 

(2.14) 
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with initial conditions 

Pi = a i ~ u l~ w l = 0 at t = 0 

and boundary conditions = 0 at z = 0, t > 0 . 

Here, 

f 1 =-[(au) x + (aw) z -<7w] 

f 2 = u p x /y ' ( u u x + w u z ) 

f 3 = a ( p z " P Y - ( u w x + w.w z ) 

= ( T P ° ' (y+U <7/2 } t - u ( p - y a ) x + ( p - y cr ) z 
+ (y ■ 1 ) (p-ycj)-(y-i)crf (x,z,t) . 


(2.15) 

(2.16) 


} (2.17 


J 


Again for convenience we shall refer equations (2.11) - (2.16) as 
problem Pp. We allow sufficient smoothness on the right hand side 
which contains terms given by (2.17) so that Pp is well-posed. In fact 
this gives us the regularity of solution of Pp. Once a numerical 
procedure is constructed for the solution of problem Pp the same 
procedure can be used to compute the solutions of Pp since 
the differential operator on both problems is the same with the bonus 
of identical initial and boundary conditions. Also once the solution (p,p^ 
), ( a, <7 j) , and ( u, ) are known, then the solutions of the nonlinear 
field are given by: 


P = * P + e 
, 2 

a - e a + e a 


(2.18) 


= e u + e _u 
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This procedure can be continued to obtain higher order 
approximations for the nonlinear problem. However, a sequence of study 
reported by Hariharan and Lester [4,5] for one dimensional problems 
and by Hariharan [6] for two dimensional problems of nonlinear 
acoustic calculations shows that only two terms are needed to 
investigate the nonlinearity, even for the case of shock waves. A natural 
question one may ask is why not solve the nonlinear problem directly, 
as in the above references, including discontinuities in the solutions 
such as shock waves?. The solutions may form shock discontinuities in 
the vicinity of the source, in which case considering two linear 
problems P^and Pjj separately will not be uniformly valid. However, we 
are interested in the sound field far away from the source, and the 
region of possible shock discontinuities is still considered as a source 
region. A full mathematical justification may be a difficult task. 


3. Formal Solutions and Estimates 


Here, we discuss the existence and uniqueness of the initial boundary 
value problems Pj and P^. We shall accomplish this by obtaining proper 
energy estimates. The first step is to write the governing equations in 
the following form: 

u. + A u + B u + C u = f (3.1) 

L X z 

where 

A = a.j, a 12 = 1. a 2 q = l/y» a^ 2 = Y ancl other a ij = 
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B = b. = 1, b^ = 1/y, b^^ - y and all other b. . = 0, 

^ J J 

C = c. . = -1, c^ = - 0^4 = i/y, C 42 = 0 the rest 0. 

J’ _ y 

Also, u= (a,u,w,p) for Pj and u = (a^,u^,v^,w^) for Pjj. 

Similarly the right hand side f has the following definition: 
f = (0,0, 0,0^ for Pj and f = (f^^jf^^)^ for P^. The boundary 
conditions are: 

w = 0 on z- 0 for Pj 

Wp 0 on z- 0 for Pjj. 


Initial conditions are : u = 0 at t = 0 for both Pj and P^. 

We want to treat the problems in the context of hyperbolic equations. Let 
us collect needed relevant information from the theory of hyperbolic 
equations. First consider the definition of hyperbolicity. Let 
and C(u) be such that 


A.(u) 


u t + 2 A i( u ) u v + C(u) u = f (x,t) . (3.2) 

j Xj 

Definition 3. 1 

If the eigenvalues of 

A (u, w) = 2 A j(u) Wj 

are real for real vectors u and w then the system (3.2) is said to be 
hyperbolic . If the eigenvalues are real and distinct, then the system is 
said to be strongly hyperbolic. 

According to this definition it is easy to verify from (3.1) that the 

2 2 A 2 2 1 

eigenvalues of Aw^+ B W£ are 0,0, (w^ + ^ 2 ' 2 anc * "^ w l + w 2 ) 2 
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satisfying hyperbolicity, but not stong hyperbolicity. The next notion we 
require is the symmetry property. Again we consider the system (3.2) 
for this purpose. In general A,(u) need not be symmetric. There are va- 

J 

rieties of procedures which are equivalent to saying the system (3.2) 

can be written in symmetric form. One of which is in the Freidrich’s 

sense; i.e., there exists a matrix valued function E(u) which depends on 

-1 

u such that the matrices B . (u) = E(u) A . (u) E(u) are 

*J J 

symmetric and then system (3.2) can be written in the symmetrc form: 


v +2B .M v = E(u) f (x,t). (3.3) 

J j 

In our considerations A and B do not depend on u implying E will not 
depend on u either. Our first goal is to obtain this matrix E which we 
call the symmetrizer. The construction follows from: 

Lemma 3.1 


There exists a matrix E such that the system (3.1) can be written in 
a symmetric form 

v +Pv +Qv - R v = E f (3.4) 

L X Z 

where, P and Q are symmetric and v = E u. 

Proof: 

The procedure consists of finding a matrix which will simultaneously 
symmetrize both A and B. The first step is to find a matrix T such that 
T *A T is diagonal. This is easily accomplished by forming the matrix T 
using eigenvectors of A. In this case T is given by 
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10 11 
0 11-1 
0 10 0 

0 0 y y 

sothatT^AT = diag (0,0, 1 ,-i) . (3.5) 

The diagonal elements are simply the eigenvalues of A . This in turn 
yields: 

0 0 0 0 
0 0 11 
0 $ 0 0 
0^00 

This matrix is not symmetric. So we further investigate possibilities of 
symmetrizing this matrix which preserves the symmetry of A. 
Consider a diagonal matrix D = diag(a,/?,y,6). This gives the above 
property for properly chosen constant diagonal elements. We observe 
that 

D -1 T _1 ATD = diag (0,0, 1,-1) , (3.7) 

which simply shows the diagonal form of A is preserved, while 

TO 0 0 0 

0 0 y/p y/P 

D *T *BTD = ® 0 0 ^ 

o Wy 0 o 

We choose <x,jS,y, and <5 in (3.8) so that the right hand side will be 





u 



symmetric. This restriction gives us the following relations: 

y/P =P/2y 

6/P - (3/2y and a arbitrary. 

Upon solving these equations, we find that one solution is P = ^2, 
Y=<*=1- Since a is arbitrary we choose it to be 1 . Then D = 
diag(i,V2,i,i). Thus the matrix G = TD gives both G ^AG and G ^BG 
as symmetric. Hence, we have the following symmetric hyperbolic 
system: 

v, + Pv + Q v - Rv = F (3.9) 
t x y 

where, v = Gu,P = GAG'!q = G A G", 1 R = - G C G' 1 and F = Gf. 

The next step is to consider the well posedness of the problems 
Pj and Pjj. Any definition of well-posedness of an initial boundary value 

problem consists of several steps. Namely, they are: 

a) Specification of spaces Hpto which F belongs, 

b) the space in which the solution v is sought, 

c) existence and uniqueness of the solution veH y for any F e Hp, 

d) continuous dependence of the solution on the function F. 

Composition of all these steps leads to a detailed analysis of the 
problem. The machinery to establish such steps follows from Friedrich 
[7], provided a suitable energy estimate is derived. Therefore, we shall 
be concerned only with deriving an estimate for these problems. It 
turns out that the energy estimate indicates the regularity of F, which 
gives a guide to modelling the source in our acoustic problem. 

To derive energy estimates for both problems, we define the 
following quantities. 
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Denote the innerproduct of two vectors by 

T 

( u,v ) = u v 

2 

Let £1 C R be the half space z ^ 0, -00 < x < 00 . Define the L 2 
norm of a vector in Q. by 

1 1 u j | ^ = / (u,u) dx dz = / / (u,u) dx dz. (3.10) 

o Q 0 -cd 

We introduce functions w = e ^ v, H = e ^ F, for some positive 
constant q. Then (3.9) becomes 

w.+Pw +Qw -Rw + fjlw = H. (3.11) 
t X z 

Now consider the derivative of the inner product 

(w,w) t = (w,w t ) + (w t ,w) 

Using equation (3.11) we have 

(w,w) t = -(Pw,w x ) - (w,Pw x ) - (w,Qw z ) - (Qw z ,w) +(w,(R+R T -2qI )w) 
+ 2(w,H ) . 

Now using the symmetry properties of P and Q derived in Lemma 3.1, 
we obtain: 

(w,w) t = - (w,Pw) x - (w,Qw) z + (w,(R+R T -2qI )w) + 2(w,H ). (3.12) 


Integrating (3.12) over Q we obtain the following energy integral: 

d_||w|| = / (w,Qw) +/ { (w,(R+R^-2qI )w)} dxdz + 

dt o -cd z=0 Q 

2/ (w,H ) dxdz . (3.13) 

Q 
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Recall our aim is to obtain an energy integral inequality. At this point 
we need the notion called ’’maximal dissipativity." A discussion of this 
concept can be found in Kreiss [8]. Suppose a boundary condition of the 
form BjW = 0 at z =0 is posed where B |is a rectangular matrix. Then 
we have the following: 

Definition 3.2: 

The boundary condition B^w = 0 at z=0 is maximally dissipative 

provided (y,Qy) ^ 0 for all y satisfying Bjy = 0. 

z = 0 1 

For the moment we shall assume there is a boundary operator B^ which 
satisfies the definition 3.2. This means that we need to prove the 
following: 

Lemma 3.2 

There is a boundary operator B^w = 0 satisfying maximal 
dissipativeness, provided w(x,0,t) = p(x,o,t) = 0. 

Remark 3.1 The resulting boundary operator B^ is exactly the 
rectangular matrix: 

0 0 10 
B l = 0 0 0 1 

Remark 3.2 In the original statement of the problem, we did not require 
p(x,0,t) = 0. This condition ensures the well-posedness of the problem 
as we will see. 

Returning to the equation (3.13) the energy integral, we use Lemma 
3.2 and obtain the following inequality, 
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d_ | |w| | £ J {(w,(R+R l - 2t\ I )w)} dxdz + 2/ (w,H ) dxdz. (3.14) 

dt o Q Q 

T 

If q Q is big enough , r) Q l > R + R , for instance, set 

1 0 = 2||R||/<S (<5<1). 

Then for q > 2 q 

-(w,(R+R T )w) + 2q (w,w) ^ 2[q - qj (w,w) ^ q (w,w) . (3.15) 


On the other hand we have the inequality, 

2 2 

(w,H) * ||w|| 1 1 H 1 1 ^ e/2 ||w|j + l/(2e) ||H||. 

Take e = q 6/2, so that 

2 2 2 
e/2 ||w|| =q6 / 4 ||w|j * q / 4 ||w|| . 

Thus the inequality (3.14) becomes 

2 2 2 2 
ilMI S- n I|w|| + n /2 ||w|| +2n'V<5 ||H|| . (3.16) 

dt o o o o 


Integrating (3.16) from time t=0 to t=T and using the zero initial 

conditions, we obtain the following inequality: 

2 T 2 T 2 

|| e v (x,T) || + q/2 f | |e ^v(x,t) 1 1 dt £ C q 1 J | |e ^F(x,t) | |dt. 

o 0 o 0 o 


(3.17) 

Here, x = (x,z), q > 2 q Q and C is a constant independent of F(x,t). 
Inequality (3.17) holds for both problems Pj and Pjj with an appropriate 
forcing function F. In summary the above procedure yields the 
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desired result. 

Theorem 3.1 

Problem Pj together with additional boundary condition p(x,0,t) = 0 is 

well-posed, i.e., for any F e there exists a unique solution v in 

L 2 (^2) * satisfying the estimate 

2 1 2 T 2 

| |v(x,T)e‘^ T 1 1 + iq f | |v(x,t)e _r)t | | ^ K q" 1 / | |F(x,t)e'' lt | | dt (3.18) 
o 0 o 0 o 

for any q > q Q >0 and for some constant K independent of F. 

Remark 3.3 Observe that u = G ^v, f = G *F, and it is readily verified 
that u satifies the same estimate as that of (3.18) with F replaced by f. 

Remark 3.4 Theorem 3. 1 suggests that the forcing function f should be 
at least in (Q) . Thus for practical considerations even if the acoustics 
is generated by pulse sources (i.e., of the "delta function" type) it should 
be approximated by a function which is in L, (&) • 

Proof of the well-posedness of the problem Pjj is similar as 
mentioned earlier and we merely state it. 

Theorem 3.2 

For any f ^ e (Q) there exists a unique solution u e l _2 (Q) such that 

2 T 2 T 2 

| |u, (x,t)e"^ T | | +!q / | |u, (x,t)e _r)t | | dt ^ K q' 1 / | |f { (x,t)e | dt 
o 0 o 0 o 

(3.19) 
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holds, provided that the additional boundary condition p^(x,o,t) = 0 is 
satisfied and the solution of Pj, i.e., u and the forcing term f are both 
in H 1 (Q). 

Remark 3.5 The requirement u,f e H ^ (Q) ( which is the Sobolev space of 
order 1) arises because forcing term for Pjj contains derivatives of the 
soution of Pjand its forcing term. 


Recall that the nonlinear solution is sought in the form: 


u = e u + e u, . 
n 1 


(3.20) 


This is as stated in equation (2.18) the linear combination of 
problems Pj and Pjj. Combining therems 3. 1 and 3.2 we have the 
following: 


Theorem 3.3 

There exists a unique nonlinear solution e (Q) for the two term 
linear solutions of Pj and P^ for the nonlinear acoustic problem 
provided that the forcing function f for Pj is in (Q). 

This theorem tells us that the smoothness of the source of acoustics 
should be more than a square integrable function. Its first derivative 
must also be square integrable. In such a situation it is necessary to 
assume sufficient smoothness on it. This becomes crucial in the 
numerical computations. If one uses a second order finite difference 
scheme, all the spatial derivatives need to be at least in (Q). Thus 
rather than considering step by step the regularity of the source, it is 
easier to approximate it by a C°° function. For example if we consider 
source term of the form f(x,z,t) = <5(x) g(z,t), where g is a smooth 
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(3.20) 


function, then one may approximate f by f where, 
f = m/Vir exp ( - m^x^). 

Similar modifications are easily made when the source is a pulse in 
the other independent variables z and t. 

To end this section let us conclude with the proof of Lemma 3.2. 

Proof (Lemma 3.2) 

Same proof holds for both Pj and Pjj. The given boundary condition in 
both problems is the normal velocity component is zero, i.e., p = 0 on 
z = 0. Therefore, let y = ( p,u,0,p) . For any vector in this form, we 
compute (w,Qw). That is 

(w,Qw) = y T G T GBG _1 Gy. 

Noting that G = TD, we have 

1 0 11 

0 V2 1 -1 

G = 

0 V2 0 0 

0 0 y y 

Then simple matrix manipulations yield 

7 

(w,Qw) = p p/y + V2 u p/y + p y. 

A variety of conditions may yield 

(w,Qw) ^ 0. However, if we choose p — 0, it is clear that 
(w,Qw) = 0 and the corresponding boundary conditions have the form: 



IB 



0 0 10 
0 0 0 1 


and the rectangular matrix is easily identified. 

Conditions such as p(x,0,t) = 0 have a physical point of view. This 
condition says that the ground is "soft." A common terminology is z =0 
is a no reflection boundary. Perhaps other boundary conditions that make 
(w,Qw) negative may correspond to other physical considerations. At 
this point we have not yet explored them. 





4. Radiation boundary conditions 

For computational purposes it is essential to truncate half space Q. 
into a finite region Q\ for example if one uses a finite difference 
scheme, then it makes computations easier if £2’ is a rectangle as 
indicated in figure 4.1. 
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Then the pieces of the boundaries F^, J^,, and F^ need to be non- 
reflecting or radiating boundaries since they must correspond to wave 
behavior at far distances. For simple wave simple wave equation such 
discussions are extensively known. A summary of these may be found in 
Hariharan [9]. In this reference, particularly the work of Engquist and 
Majda [10] is noted. What follows is an attempt to extend the idea in 
[10] to obtain necessary boundary conditions for the linear problems 
under consideration. For this purpose we shall be concerned with only 
the problem P^ . The same radiation conditions are applicable to 
problem Pjj. Recall that problem P^is prescribed by equations (2.5) - ^ 
(2.10) together with the additional boundary condition p(x,0,t) = 0, at 
z = 0. Suppose we are interested in the radiation boundary condition 
on the boundary I^.It is sufficient to treat this boundary alone 
for obtaining radiation conditions. That is to say treat the problem as a 
half space problem, with F^ playing the role of of the x axis ( -oo < x < 
oo). 


To follow the idea of [10], we take Laplace transform with respect to 
time t (since the initial values are specified ) and Fourier transform 
with respect to x of equations (2.5) - (2.8). Using initial conditions 
(2.9) we obtain. 


A A A A 

sa + i(u+w z -w = 0 

A A 

s u + i£/y p = o 

A A 

s w + 1/y p z - = 0 

A AAA 

s p + iy£ u + y w - w =0. 


(4.1) 
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At far distances the effect of the forcing term f vanishes. Thus we seek 
homogeneous solutions of (4.1) with dependence exp(Az). This leads to 
the characteristic equation 

s i£ A-i 0 

0 s 0 i£/y 

det = 0. (4.2) 

1/y 0 s (\-i)/y 

0 iy£ (yA-1) s 

This has roots 

X 1>2 = *±,i( s,£) # (4.3) 

where 

Ms,?) = rt + s 2 + ? 2 + f, 2 ! 2 /* 2 ) i (4-4) 

with 

p 2 = (y 2 -t)/r- ( > 0) (4.5) 

For a decaying wave we choose the negative root of (4.3). Moreover, 
jj{ s,£) can be written in the form 

Ms,?) = (s 2 +w 2 (?)) 5 (s 2 +w|(?))Vs (4.6) 

Indeed one can obtain a theoretical solution of (4.1) by variation of 
parameters using the homogeneous solution dictated by exp(Xz) with the 
values of X given by (4.3). But the dificulty will be to invert the 
transforms using the boundary conditions. A similar difficulty arises 
at far distances even without the source term. It is easy to see 
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that the integrand will contain terms of the form exp(^z - jj( s,£)z ). 
From equation (4.6) we obtain and explicitly as follows: 

0J[K) = il(i+? 2 +2«) J + u H 2 -2^)h 

« 2 (f) = MHe 2 +2/S?>* + (1 +f 2 -2 wh . 

Let s = ir. Then (4.6) becomes 

|u(ir,f) = l/(ir) ((r 2 -« 2 (Q) (r 2 -w^ © ) 1 ^ . 

For propagating waves we require that /j to be imaginary. This is 
ensured by |r| > Wj(£) and |r| < From (4.7) one finds such a 

requirement is satisfied provided |£/r| << 1 * or equivalently |£/s| <<1. 
Thus we shall be concerned with approximating /4s, £) for large values 
of s. To see this let us emphasize that the solution of X we seek is of 
the form 

A = i - /4s,£). (4.9) 

Multiplying equation (4.9) by exp(Az), we obtain an associated 
differential operator 

= i p - ,i(s,f) p . (4.10) 

In the pseudo differential operator terminolgy /4 s,£) is the symbol of an 
associated pseudo differential operator.. To obtain radiating solutions 
equation (4.10) needs to be inverted for both the Laplace transform and 
the Fourier transform. Indeed a perfectly absorbing boundary condition 
arising from (4.10) is given by the inversion, 


(4.7) 


(4.8) 
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3 100 00 ^ . . £. iuo (JU ^ t 

3 z/ / pCz;^,s) e s ^ x d£ds = / / p(z;£,s){$ - /4s,£)} e st ^ x d£dx. 


dz J . 

-loo -oo 


100 00 ^ 
r / f 

-ioo -oo 

The above expression simplifies to 
ioo oo 

r / / 

-ioo -co 


3 ioo oo ^ . ... 

gg = 4 P - / / p(z;?,s) e st ‘ * x d£ds . (4.11) 


Equation (4.1 1) provides a boundary condition that is nonlocal in time 
and space. A similar procedure to that discussed above will hold for 
boundary conditions on the boundaries T ^ and ^ provided we take the 
Fourier transform with respect to z and construct differential operators 
in the direction of x. Boundary condition (4.11) is not easy to 
implement. However, if we approximate the symbol /i(s,£) for large 
values of s then it is possible to obtain approximate local boundary 
conditions from (4.11). To do this we consider /i(s,£) again and 
investigate its nature when |£/s| << 1. We rewrite fj as 

/i(s,0 = s[ 1 + a+| 2 )/s 2 + jS 2 f 2 /s 4 ] 4 . (4.12) 

A crude approximation is /j - s. Substitution of this approximation in 
equation (4.1 1) gives the boundary operator 

(4,13) 


This is a possible boundary condition. Using the Taylor approximation of 
(1+x)^ for small x, we see from (4. 12) that the next level of 
approximation is 




23 



H.14) 


,4s, fl = s[ 1 + l/(8s 2 ) + i ? 2 /s 2 ] . 

We differentiate (4. 1 1) with respect to t and substitute the 
approximation (4.14) to obtain the next order boundary condition 

P zt = J P - (P lt + P/8 - i P XX ) ' (4 - 1S) 

Similarly, higher order accurate boundary conditions can be derived. 
This process seems elegant. However, not all such boundary conditions 
yield stable results. That is to say, well-posedness of the problems is 
not guaranteed with the ail such boundary conditions. At this point it 
remains to be shown that we can derive energy estimates of the form 
(3.18) with boundary conditions of the above type on the boundaries 
IT,, and [^.Discussion of these results including corresponding discrete 
versions of our problems will be reported elsewhere. 
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